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Unitary transformations are an essential tool for the theoretical understanding of many systems 
by mapping them to simpler effective models. A systematically controlled variant to perform such 
a mapping are perturbative continuous unitary transformations (pCUT) among others. So far, this 
approach required an equidistant unperturbed spectrum. Here we pursue two goals: First, we extend 
its applicability to non-equidistant spectra with the particular focus on an efficient derivation of the 
differential flow equations, which define the enhanced perturbative continuous unitary transforma- 
tion (epCUT). Second, we show that the numerical integration of the flow equations yields a robust 
scheme to extract data from the epCUT. The method is illustrated by the perturbation of two- leg 
spin ladders around the strong-rung-coupling limit for uniform and alternating rung couplings. The 
latter case provides an example of perturbation around a non-equidistant spectrum. 
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I. INTRODUCTION 



from the Hamiltonian 



Quantum many-body systems with correlations are no- 
toriously difficult to describe theoretically. Many an- 
alytical and numerical tools have been developed to 
tackle such problems. A tool which is employed ubiq- 
uitously are unitary transformations. Famous applica- 
tions are the fermionic Bogoliubov transformations in 
the mean-field theory of superconductivity by Bardeen, 
Cooper and Schrieffer^ or the bosonic Bogoliubov trans- 
formations arising in linear spin wave theory of quan- 
tum antiferromagnets 2 . These are exact transforma- 
tions which use the algebraic properties of fcrmions and 
bosons, respectively. They yield diagonal Hamiltonians 
if they are applied to bilinear initial Hamiltonians. 

Another class of unitary transformations are those 
which are not exact but approximate because they rely 
on an expansion in a small parameter. A well-known 
example is the antiferromagnetic Heisenberg exchange 
coupling J as it is derived from a half-filled Hubbard 
model with hopping t and local repulsion U implying 
J = 4t 2 /U, see for instance Ref. |31 Obviously, higher 
contributions O (t 3 /U 2 ) are neglected. But they can also 
be computed systematicallj^^. 

Moreover, the Hubbard model is not diagonalized 
by the transformation but mapped to an effective spin 
model. This mapping implies a simplification because the 
relevant part of the Hilbert space (here: spin degrees of 
freedom) has been separated from the remainder (charge 
degrees of freedom). The remainder does not need to be 
considered anymore. It is said that it has been eliminated 
or integrated out. 

Another famous example in the same line is the 
Frohlich transformation^ which eliminates phononic de- 
grees of freedom from an electron-phonon system in lead- 
ing order of the coupling to derive an electron-electron 
interaction from an electron-phonon coupling. Starting 
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this transformation generates an attractive interaction in 
the BCS-channel 
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where q := k — k, Ae = This explains the forma- 

tion of Cooper pairs and conventional superconductivity. 
It is interesting to note that in standard treatments the 
interaction is usually approximated by a constant, leav- 
ing out any discussion of the resonance singularity in (|3| . 

It is, however, possible to achieve the elimination of 
the phonon degrees of freedom by a different, continuous 
unitary transformation (CUT). This approach relies on a 
continously parametrized antihermitean generator rj(£) = 
—rf{£) of the differential unitary transformation 



d e H(£) = [ v (£),H(£)} 
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of the Hamiltonian H(£); the transformation starts at 
£ = and endfPHIH a t £ = oo. 

One possible choice for the generator leading to a con- 



vergent flow^ for £ — > oo is r\ 
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[Hd,H] where 7?d 



is the diagonal part of the Hamiltonian. Integrating the 
flow equation Q from £ — to oo yields for the BCS 
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The eye-catching fact in V^p is that it does not have a 
resonant energy denominator. Hence this result is much 
smoother. In particular it implies an attractive interac- 
tion for all parameters. 

The standard BCS interaction is a constant up to some 
phononic cutoff energy WDcbyc- This result can be derived 
rigorously by a modification of the generator. In an eigen 
basis of Hd the matrix elements of ?y MKU are choserP^nSS 
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where Q(x) is the Heaviside step function. Again, there 
is only attractive interaction. In addition, the interac- 
tion is only active in a restricted energy interval and zero 
outside. 

It is very remarkable that all three approaches [3j [5] 
and [6] are different in their outcome although they do 
the same: Eliminating the linear electron-phonon cou- 
pling. We stress that this is not a spurious result, but 
relies on the fact that the unitary transformations are in- 
deed different even in leading order. They express virtual 
processes in a different way. But the energy conserving 
processes at Ae = are the same in all three results. 
This has to be so because such scattering processes can 
in principle be measured which implies that they have to 
be independent from the chosen basis. 

We include the above results to demonstrate the power 
of the freedom to choose an appropriate basis. The 
derivation of the BCS electron-electron interaction illus- 
trates the versatility of continuous unitary transforma- 
tions. Moreover, it exemplifies that even a seemingly 
simple task, the derivation of an effective interaction in 
leading order, provides unexpected surprises. We advo- 
cate to do it in a less singular way with effective interac- 
tions which are restricted in energy range. 

The main goals of the present article are twofold. First, 
we show how CUTs can be used to perturbatively de- 
rive effective Hamilton operators in real space. This goal 
has been realized for an unperturbed Hamiltonian with 
equidistant spectrum by perturbative CUT (pCUT^pHio]. 
The gist of the pCUT is recalled below. In the present 
work we enhance the applicability of such an approach 
to unperturbed non-equidistant spectra by formulating 
the CUT directly in second quantization. The resulting 
transformation will be called enhanced perturbative CUT 
(epCUT) for distinction. The approach is exemplified for 
a uniform and for an alternating spin ladder. The latter 
has a non-equidistant spectrum if only the rung couplings 
are considered. 

The second main goal is to establish a robust extrapo- 
lation of the perturbative results of the epCUT. We will 



show that a direct evaluation of the perturbatively es- 
tablished flow equation provides a very robust and reli- 
able way to extrapolate the perturbative results. This 
approach will be called directly evaluated enhanced per- 
turbative CUT (deepCUT). 

The article is set up as follows. In the remainder of 
this section we briefly introduce the perturbative CUT 
and the self-similar CUT (sCUT) as predecessors of the 
epCUT and the deepCUT. In Sec. II, we introduce our 
paradigm model, spin ladders, for which we illustrate the 
general approaches. In Sec. Ill, we derive the epCUT 
and develop the deepCUT from it. Many technical as- 
pects are discussed; a focus are simplifying rules which 
allow us to compute high orders efficiently. In Sec. IV, 
results of the epCUT and the deepCUT are presented for 
the uniform antiferromagnetic spin ladder with S = 1/2. 
Results for the alternating spin ladder, which does not 
have an equidistant unperturbed spectrum, are shown in 
Sec. V. The article terminates by the conclusions in Sec. 
VI. 



Perturbative Continuous Unitary 
Transfer mat ion 



We draw the readers' attention to the fact that we are 
dealing from now on with CUTs with a unique reference 
state. This means that the ground state is mapped by the 
CUT to the vacuum of excitations. In the Introduction, 
the mapping to effective models such as the Heisenberg 
exchange model or the BCS-model still left a many-body 
problem to be solved. 

High order series expansions have long been used to 
compute reliable ground state energies and dispersions 
in strongly-correlated systems^. These quantities can 
be computed straightforwardly because the states are 
uniquely determined by their quantum number, for in- 
stance the momentum, even before the perturbation is 
switched on. This means it is sufficient to perform 
the perturbation for a one-dimensional subspace of the 
Hilbert space. States of two and more particles are more 
subtle because their subspaces are extensively large. For 
instance, binding energies cannot be computed as series 
unless the binding occurs already in linear order. Gener- 
ally, unitary or orthogonal transformations must be in- 
troduced to defi ne the perturbative approach on large 
subspacePESEDESI. 

The perturb ative CUT was the first approach to 
achieve this goa l 18 * 19 *. Its starting point is a Hamiltonian 
which can be written in the form 



H = H 
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where Hq is the unperturbed Hamiltonian with an 
equidistant spectrum. For simplicity we set its energy 
spacing to unity. Each energy quantum can be seen 
as an elementary excitation, a quasi-particle, so that 
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H = counts the number of quasi-particle up to an ir- 
relevant offset. The expansion parameter is x and the 
terms in the perturbation are split according to their 
effect on the quasi-particle number Hq: The terms in 
T m increase the number of energy quanta by m. Obvi- 
ously, T_ m = Tj n holds. Generically, there is an upper 
bound N > \m\ to the change of energy quanta. Choos- 
ing 77 MKU = J2m=-N s E n { m )T m a flow is defined via Q 
which generates a hierarchy of differential equations in 
powers of a;. In each finite order, the differential equa- 
tions are closed and can be solved by computer aided 
analytics. Eventually, one obtains the general expansion 
for an effective Hamiltonian which conserves the number 
of elementary excitations 

oo 

H cB = H + J2 £ C{rh)T(m). (8) 

k=l dim(m)=fc,M(m)=0 

The components of the vector fh take the values 
—N, -N + 1, . . . , N - 1, N; the vector has the dimen- 
sion dim(m) = k. The coefficients C(m) are sim- 
ple real fractions and T(m) stands for the product 
T mi T m2 . . . T mkl T„ lk 19 . The conservation of the number 
of quasi-particles is implied by the cross sum M(fh) = 
where M{m) = Sj=i 171 j- 

The result ^ is very general; to put it to practical use 
its irreducible effect on zero, one, two, and more quasi- 
particles is computed. In this way, the effective Hamilto- 
nian is obtained in second quantized forrrpSI. Remarkable 
achievements of this approach are a quantitative under- 
standing of inelastic scattering in spin ladders^^U, of 
spectral dens ities in spin chains^!, of excitations in the 
Kitaev modePES f excitations in the toric code^EU, 
and of the ionic Hubbard modeP^to name a few extended 
systems where the ground state is described as a vacuum 
of excitations. 

Conceptually, the most significant achievement of 
pCUT is that whole subspaces are treated perturbatively. 
The generality of the pCUT result ^ is surely one of its 
advantages. The fact that it can only deal with equidis- 
tant spectra is a certain caveat. Another caveat is that 
the approach does not allow for modifications of the gen- 
erator. 



B. Self-Similar Continuous Unitary 
Transformations 

One way to circumvent the above mentioned restric- 
tions concerning the unperturbed spectrum and the 
choice of the generator is to pass from a perturbative 
evaluation to a self-similar one. The approach follows a 
straightforward strategy. One chooses a set of operators 
which serves as a basis. The Hamiltonian and the gen- 
erator are described as linear combinations of these op- 
erators. By commuting Hamiltonian and generator and 
re-expanding the result in the same operator basis the 
flow equation (HI) induces a differential equation system 
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Figure 1. (Color online) Schematic representation of the sym- 
metric (alternating) S = V 2 Heisenberg ladder in the thermo- 
dynamic limit. 

(DES) in the coefficients of the basis operators. A more 
detailed description follows in Sec. Ill below. We stress 
that in the latter step a certain truncation is required^. 
Unless the set of operators is closed under commutation, 
the commutator [77, H] comprises terms which cannot be 
expanded exactly in the operator basis. Thus this step 
generically requires an approximation. The Hamiltonian 
is kept in a self-similar form defined by the selected op- 
erator basis. 

Depending on the system the truncatin g appro xima- 
tion can be controlled by a small paramete r 15 * 38 * 39 * or by 
the spatial locality of the selected set of operator d 7 * 40 [ 
We stress that in the sCUT approach the choice of the 
operator basis and of the generator uniquely defines the 
DES of the flow equation. Clearly, the advantage of the 
sCUT over the pCUT is its larger versatility. Yet it is 
less general in the sense that the flow equation has to be 
solved for each model and each operator basis anew. 

To derive a systematic perturbative expansion by 
sCUT is not an obvious step. It is one of our two main 
goals to show how this can be done and how it can be 
done efficiently. Thus the derivations and considerations 
in Sec. Ill are based on the sCUT approach and combine 
it with a perturbative expansion in order to reach the 
enhanced perturbative CUT. 

II. MODEL 

To illustrate the performance of the (de)epCUT we 
consider the S=Y2 antiferromagnetic two-leg Heisenberg 
ladder (uniform spin ladder) and an extension with an al- 
ternating rung coupling (alternating spin ladder) as test- 
ing ground, see Fig. [l] The Hamiltonian reads 

H = JJ_H1 + Jj_H1 + J\\H\\ (9a) 
L/2-1 

Hl= Y, S 2r- S l (9b) 

r=0 
L/2-1 

Hl= £ S L 2r+1 -Sl +1 (9c) 

r=0 
L-1 

=J2( S r- S^+i + S? • S« +1 ) , (9d) 

r=0 
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where r£Z. The rung number is denoted by r and the 
legs by L and R. We define the ratio between the leg 
coupling J ii and the even rung coupling JJ_ as x := 
and the ratio between the odd rung coupling and the 
even rung coupling JJ_ as y := 

In the limit of JJ_ = J]_, i.e., y = 1, the Hamilto- 
nian describes the uniform spin ladder. This model has 
been subject of intensive studies, see Refs. H7J Iffl and 
|42] and references therein. Thus it constitutes a suitable 
reference model to test the epCUT. It has been investi- 
gated by several di fferen t methods, such as density ma- 
trix renormaliza tiorP^ ^, exact diagonalizatiorPH, contin- 
uum field theorjifiEa qua ntum Monte Carlc@^, high or- 
der series expansion s'^, in cluding methods ba sed on 
CUTs, such as sCU T 4 " ! 51 ' 52 ' and pCUT™5EZE3] If the 
results of the epCUT agree with this data, the efficiency 
of the epCUT for the expansion around an unperturbed 
equidistant spectrum is verified. 

To illustrate that the epCUT represents an advance- 
ment compared to pCUT we will show results for the al- 
ternating spin ladder as well. This system does not have 
an equidistant spectrum because the rung couplings are 
not equal Jj_ ^ Jj_. Hence it cannot be dealt with by 
pCUT. Without loss of generality we consider J° > Jj_ 
implying y > 1. 

For the alternating spin ladder we expect a lowering 
of the ground state energy upon rising y because the 
expectation value of (S r • S^) is negative. The unit cell 
includes two rungs which implies two triplon branches in 
the Brillouin zone (BZ). For y = 1 the branches meet at 
the BZ boundary (k = i^/a). For y > 1 a band gap of 
the order of \y — 1| opens at k = ±*/2 separating the two 
band. 

To define a starti ng p oint for the CUT the bond oper- 
ator representation^ 4 ^ [ s used. A possible eigen basis of 
the local operators S r , is given by the singlet state 



hardcore-boson commutation relation 

7 

(15) 

The el ement ary magnetic excitations (S=l), known as 
triploniPHHl can be continuously linked to the local 
triplets. 

Represented in second quantization in terms of the 
triplon creation and annihilation operators the Hamil- 
tonian reads 



H 
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where 



H ± = -7 E 1+ E f lrt a , 



r—2a r—2a 



(17a) 

H°X = -- A E 1+ E *a,r*o,r (17b) 
r=2a+l r=2a+l 

^11 = 2 E/ i^a.r^a,r+l + ^a,r+l^cx,r J (l? c ) 

+ - V J t f t t 
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where a, r € Z. This form of the Hamiltonian enters all 
the calculation described in the following. 
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III. DERIVATION 
A. Flow equation in second quantization 



and the three triplet states 

4 \s) ■= K) = ^ (|tt) - III)) (ii) 

4| S ):=|^> = ^=(|tt) + l44}) (12) 
tt\s):=\t z ) = ^=(\U) + \m- (13) 
For x = 0, the ground state of the system is given by 

|0>:=ni s )r- (14) 

r 

This vacuum of triplets serves as our reference state. The 
local operators Q, r , t* and t\ r (t xr , t y r and t z r ) create 
(annihilate) an excitation on rung r. They satisfy the 



Simil ar to th e implementation of previous CUT 
methods ^ 12 * 20 * 40 !, we formulate the flow equation Q for 
the coefficients of the monomials {A{\ of operators in 
second quantization. The Hamiltonian is parametrized 
by 

H(t) = Y / h l (£)A l (18) 

i 

with the ^-dependent coefficients hi(£). The generator 
reads 

m = E ^ : = E wwim (19) 

i i 

with f] being a superoperator denoting the application of 
a particular generator scheme such as those discussed in 



Ref. HOI Expanded in the operator basis {A{\ the flow 
equation Q reads 

Y dMQAt = Y hj(£)h k (£) [r)[Aj],A k ] . (20) 

i jk 

Comparing the coefficients of different monomials, the 
flow equation Q becomes equivalent to a set of ordinary 
differential equations for the coefficients hi(£) 



d e hi(£) = Y D ijkhj(£)h k {£). 



(21a) 
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The commutator relations between the basis operators 
are encoded in the coefficients D^k of the bilinear differ- 
ential equation system (DES). These coefficients Dij k are 
in general complex numbers. For the spin ladders under 
study they are given by integers or fractional numbers. 
We call a single Dijk a 'contribution' of the DES. The 
contributions are obtained from 



Y D ikjAi = {r)[Aj],A k } 



(21b) 



by comparing the coefficients of the expansion of the com- 
mutator monomial by monomial. 

In this way, the problem of solving the flow equation 
is transformed into the algebraic problem of calculating 
the coefficients of the DES (21b) and of the subsequent 
numerical solution of Eq. 



(21a 




Figure 2. Sketch of the epCUT algorithm to calculate the DES 
for the iterative calculation of diH^ A \ Due to the commu- 
tators [rj W ,H (3) ], . . . , [?7 (3) ,# (1) ], new terms with O min = 4 
emerge. Thus, the calculation of the block [r]^, H^] has to 
be carried out at last and self-consistently because it generates 
monomials contributing to the generator in the same order. 
If H is not (block-)diagonal, both [r/ 4) , H {0) ] and [r/ m , # (4) ] 
have to be calculated simultaneously in a single self-consistent 
loop. 



B. Perturbative expansion of the flow equation 

Here we consider the perturbative solution of the flow 
equation which yields the resulting effective Hamiltonian 
in the form of a perturbative series. Henc e this solution 
generalizes the established pCUT approac h 1 18 * 19 !. To this 
end, we decompose the initial Hamiltonian 



H = Hn 



xV. 



(22) 



into an unperturbed part Hq and a perturbation V. In 
contrast to pCUT^, we do not require the unperturbed 
part to have an equidistant spectrum. We do not require 
other restrictions either although we will see that the 
method works best for a (block-)diagonal Hq. 

We aim at the perturbation series up to and including 
order n in x. Thus we expand the flowing Hamiltonian 



H{£) = Y ff(m) > H 



(m) 



m=0 



into terms of order x m up to m < n. Expanding the 

in the operator basis {A{\ we perform the expansion in 

powers of x by expanding the coefficient hi(£) of Ai 
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At / = 0, the initial values fi m \o) are fixed by the initial 
Hamiltonian ( 22 ) and its representation in terms of the 
{Ai}. Applying ( 24 1 to Eq. (21a) one obtains 



di Y ^ 



(m) , 



YD ijk y * p+ vj p) w/i 9) w. 

j,k p,q=0 

(25) 



For the prefactors of x m this implies 

j,k p+q=m 



(26) 



We stress that the contributions Dij k do not depend on 
the order m of the coefficients, but only on the algebraic 
relations between the corresponding monomials. Hence 
they need to be calculated only once. Moreover, Eq. 



(23) (261 defines a hierarchy between the coefficients because 



f- {£) is influenced only by coefficients of the same order 
m or lower, but not by coefficients of higher orders. 



C. Algorithm 

A key task in the implementation of epCUT is the de- 
sign of an efficient algorithm to identify the monomials 



G 



and to calculate exactly the commutators which are rel- 
evant for the transformed Hamiltonian in the order of 
interest n. Henceforth, we call the order we are aiming 
at the 'targeted' order. 



Based on Eq. ( 26 ) we can calculate each order m based 
on the results of lower orders. Order zero is trivially given 
by the representation (22 1 if t][Hq\ = which means that 
H is block-diagonal. The calculation of the commuta- 



(m-l) 



can be carried out 



According to Eq. (21b I, the 



tors [^\H { ' 
independently, see Fig 

commutator [^[A,-], A^] can be written as linear combi- 
nation of monomials of which the prefactors define 
the contributions -Djjfc of the DES. For those monomials 
not yet present in the Hamiltonian, a new monomial has 
to be included in the operator basis with a unique index. 
We call the order in which a monomial occurs for the first 
time its minimum order O min (A;). 

We stress that in the evaluation of i/' 9 '] , the 

commutator [^[Aj], Aj,] needs to be calculated only if 
O nlin (Ai) = p and O min (Aj) = q. For all monomials 
with lower O min (A) and/or lower O min {Aj), the com- 
mutators have already been calculated in lower orders. 

The calculation of the commutators for [r]^ 71 ', ifW] is 
special because its result may include new monomials of 
the same minimum order m which were not considered 
so far. Since these monomials also enter the commuta- 
tor via rf m \ the block [7/ m ),iJ( )] has to be iterated 
until no new monomials occur: Then self-consistency 
is reached. This should be done once the inner blocks 
[^p >0 \ ij(«>°)] are finished. 

If the unperturbed Hamiltonian H is local, the com- 
mutation of monomials from ?y( m ) and Hq lead to mono- 
mials acting on the same local cluster or smaller subclus- 
ters. Furthermore, if the local Hubert space of the cluster 
is finite, the number of new monomials which can be gen- 
erated by iterative commutations with Hq is bounded by 
the finite number of linearly independent matrices on this 
finite-dimensional Hilbert space. Then the iterative loop 
is guaranteed to terminate after a finite number of cycles. 
In the symmetric ladder model (see Eq. (9dl) the local 



Hilbert spaces are finite so that a finite number of cycles 
is certainly enough. Even better, the commutation of the 
monomials in terms of triplon creation and annihilation 
operators with Hq = J^H^ + does not generate 

any new monomials so that no iterations are needed in 
the calculation of [n^ >0 \ #(9> )] . 

If the unperturbed Hamiltonian Hq has also non- 
(block-)diagonal terms, the generator includes terms of 
order zero. Therefore, the blocks H^] have to 

be evaluated self-consistently as well. Since any term of 
the Hamiltonian may also appear in the generator, the 
blocks [rf m \H^] and [n^ Q \H^] have to be calculated 
simultaneously within a joint self-consistency loop. 

For the sake of completeness, we note that in the spe- 
cial case Hq := H, i.e., considering the total Hamilto- 
nian as the unperturbed one, the whole algorithm con- 
structing the DES reduces to the calculation of the block 
W ), Jj(°)l . This has to be done self-consistently with 



respect to both the generator and the Hamiltonian. This 
approach is the o ne em ployed in the self-similar CUT 
(sCUT) previousljl29M0|. Since for H = H "the un- 
perturbed" part Hq includes non-local and non- (block- 
diagonal terms and perhaps refers even to an infinite 
local Hilbert space for any but the simplest models, the 
iteration of commutators will not terminate. Thus addi- 
tional truncation criteria are needed whose validity needs 
to be justified. 



D. Perturbative evaluation 

To evaluate the perturbation series for the ground state 
energy or the dispersion relation of a concrete system, the 
first step is to write the Hamiltonian in second quanti- 
zation and to identify the relevant monomials. We illus- 
trate the approach by the second order calculation for 



the symmetric spin ladder in Eq. (9d) using the quasi- 
particle conserving generato j 19 * 40 l 

The operator basis {A{\ is given in Tab. [I] with Aq and 
Ax for the terms in Hq (O m i n = 0) and A 2 to A 5 for the 
terms in V (O m in = 1)- We combined certain monomi- 
als whose prefactors must be the same due to symmetry 
and/or hermiticity into one element of the operator basis 
Ai (cf. Sec. HIE 2 1. The advantage is that less operators 



need to be tracked. The algorithm is not affected by this 
step except that the comparison of coefficients is a bit 
more complex. 

Following the algorithm described above, the commu- 
tators of the block [77W, HW\ are calculated to complete 
the first order. The contributions to the DES obtained by 
comparison of coefficients are given in Tab. [H] Then, the 
contributions in second order are evaluated in the blocks 
[7/ 1 ), ii^ 1 )] and [r]^ 2 \H^] leading to the new basis op- 
erators Ag_i7 with O m i n = 2. 



Next, the perturbative flow equation (26) has to be 
solved. We do this numerically using a standard fourth 
order Runge-Kutta method^. The initial values for the 
coefficients in different orders of x are read off the initial 
Hamiltonian. They are zero for all basis operators and all 
orders which are not present in the initial Hamiltonian. 
We use a basis of only normal-ordered operators except 
for Aq = J2 r ^ so that the series expansion of the ground 
state energy per rung Eq is obtained in the limit of £ — > oo 
from the prefactor of Aq 



Eo=J2 f™(oc)x m + O {x n+1 ) (27a) 

m— 

3 3 



(27b) 



Note that this result requires only three equations in the 
DES. 

Likewise, the dispersion relation is determined from 
the renormalized coefficients of the hopping terms Ai,A§ 
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Table I. Basis operators Ai (simple combinations of monomials obeying symmetry and/or hermiticity) occuring in a second 
order epCUT calculation for the uniform spin ladder using the particle conserving generator scheme rj pc - The third column 
contains the initial coefficients hi(£ = 0), the fourth the final renormalized coefficients hi{t = oo). The minimum order O min 
is given in which the corresponding operator occurs for the first time; OSS^ is the highest relevant order of the coefficient for 
computing the ground state energy and is the highest relevant order for computing the dispersion. The terms marked in 

light gray are irrelevant for the computation of the ground state energy in second order; the terms in dark gray are irrelevant 
if the dispersion is computed. 
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Table II. Non-vanishing contributions Dijk to the differential equation system (DES) of the uniform spin ladder using the 
particle conserving generator scheme. Operators and contributions marked in light gray are irrelevant for the computation of 
the ground state energy in second order; those in dark gray are irrelevant for the dispersion. 



and As 

n 

uj{k) = (A m (oo)x m + 2/ 5 m (oo)a; m cos(fc) (28a) 

m=0 

+ 2/ 8 m (oo)x m cos(2fc)) + O (x n+1 ) (28b) 

3 1 

= 1 + -x 2 + x cos(fc) - -x 2 cos(2fc) + O (x 3 ) , 

(28c) 

which require five equations, only two more than the 



ground state energy. 



E. Optimizations 

The epCUT method presented so far can be applied to 
a wide range of models in order to calculate a perturba- 
tive expansion of decoupled quasi-particle spaces. With 
increasing order the number of representatives in the ef- 
fective Hamiltonian, the runtime and the memory con- 
sumption rise exponentially, see Fig. [3] One is interested 




2 4 6 8 10 12 14 16 6 8 10 12 14 16 

order order 



Figure 3. Panel (a): Number of representatives in the effective Hamiltonian vs. the order of the calculation for various 
optimizations aiming at the ground state energy using 770 and all symmetries. Highest to lowest curve: full Hamiltonian, 
basic simplification rule, extended rule, full reduction of the DES based on the exact O max . Panel (b): Runtime time for 
the construction of the DES vs. the order of the calculation with more and more optimizations using rjo and all symmetries. 
Highest to lowest curve : full Hamiltonian without simplification, basic a-posteriori simplification rule, extended a-posteriori 
rule, additional use of the basic a-priori rule, additional use of the extended a-priori rule. The computations were done on an 
Intel Xeon CPU (E5345, 2.33GHz, single thread). 



generator order # representatives runtime RAM 

scheme [dd:hh:mm] [GB] 

On 17 51,731,694 2:17:14 8~X 

0:n, l:n 15 107,513,297 13:09:12 17.3 

0:n, l:n, 2:n 13 51,371,642 11:09:47 8.0 

Table III. Number of representatives in the operator basis, 
total runtime and memory consumption for various generator 
schemes in the highest order calculated. The computations 
were done on an Intel Xeon CPU (E5345, 2.33GHz, single 
thread) with full optimizations. 

in increasing the order of the calculation as high as possi- 
ble because this generically enhances the accuracy of the 
calculation: More and more orders kept imply that more 
and more physical processes with an increasing range are 
taken into account. 

To increase the order, more efficient generator schemes 
and the symmetries known from sCUT can be exploited. 
Focussing on selected quantities of interest, the pertur- 
bative foundation of epCUT allows us to optimize the 
algorithm even further. Generic performance data possi- 
ble with full optimizations are given in Tab. |III| 

In practice, every optimization is carefully checked by 
comparing the results of the optimized faster program to 
the results from the slower program before optimization. 
In this way, one can be sure that no errors are introduced 
by incorrect assumptions. 

1. Generator scheme 

The (quasi-)particle-conserving generator scheme fj pc 
used in our example decouples all subspaces of differ- 



ing numbers of excitations, i.e., quasi-part icles, and sorts 
them in ascending order of their energ^E H 19 ! 40 !. In most 
applications, however, only the ground state and the low- 
lying excitations are of interest. Consequently, the com- 
putational effort can be reduced by choosing a more ef- 
ficient generator scheme which targets the quantities of 
interest only. In 2010, Fischer, Duffe and Uhrig^ pro- 
posed a family of generator schemes based on modifica- 
tions of rfp C where only the first q quasi-particle spaces 
are decoupled from the Hilbert space. The corresponding 
generator reads 

q 00 

rj q [H(£)] := £ £ (h){1) - h{(£)) . (29) 

j=0 i=j+l 

In this notation, ifj comprises all monomials of the 
Hamiltonian creating i and annihilating j quasi-particles. 
For instance, the ground state generator 

rfo[H(e)}=J2(Hi(i)-HUe)) (30) 

i 

incorporates monomials which consist purely of either 
creation or annihilation operators. Compared to the full 
quasi-particle conserving generator, the effort to compute 
the corresponding DES is reduced significantly. In anal- 
ogy to ?7p C , the decoupled quasi-particle spaces are sorted 
according to energy. Thus the ground state energy is 
given by the vacuum energy of the effective Hamiltonian 
H(oo). If additionally the dispersion is calculated, the 
one-quasi-particle subspace has to be decoupled using 
rji. For decoupling higher quasi-particle spaces, analo- 
gous generator schemes can be used. But the increase in 
efficiency compared to the full quasi-particle conserving 
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generator becomes less and less significant because the 
generator schemes rj q do not conserve the block-band- 
diagonal structure of the Hamiltonian in contrast to rj pc . 



2. Symmetries 

For models defined on infinite lattices, it is necessary 
to use the translation symmetry in order to be able to 
work directly in the thermodynamic limit. In addition, 
the presence of other symmetries leads to linear depen- 
dences of coefficients of monomials which are linked by 
the s ymm etry transformations of the Hamiltonian. As in 
sCUTEES] this redundance can be significantly reduced 
by passing from simple monomials to symmetric linear 
combinations of them. Each of these polynomials is in- 
variant under symmetry transformations of the Hamil- 
tonian and requires only one prefactor where the single 
monomials would need many more. In our example (Tab. 
|l|) and in the following calculations the size of the oper- 
ator basis is reduced by a factor of almost 24 exploiting 
self-adjointness, reflection and spin symmetry. 



3. Reduction of the differential equation system (DES) 

Targeting only certain quantities up to order n, such 
as the ground state energy or the one-particle dispersion, 
the DES can be reduced. Here we discuss how this can 
be done in practice. 



Besides the minimum order O n 



a maximum order 



O max can be a assigned to each monomial and its coeffi- 
cient hi. The maximum order is the highest order of the 
series of hi which still has an influence on the targeted 
coefficients up to order n. For instance, complicated pro- 
cesses involving many quasi-particles do not influence the 
ground state energy directly, but only via other process. 
Then their O max is much lower than the targeted order 
n. Technically, this is due to the hierarchy of the DES 



( 26 ) which implies 



O max (^) > O max (A) - O min (A k ) 

c(^i) — O min (Aj), 



(31a) 
(31b) 



where the equality holds if we consider only a single con- 
tribution Dijk 7^ 0. The inequality takes into account 
that there may be many pairs (i, k) for a given j. Thus, 
O max (Aj) is the maximum value of all those right hand 
sides 



Omax(Aj ) 



max [O n 

{i,k\D ijk ^0} 



(Ai) - O min (A k )} . (32) 



If Aj is targeted, for instance the ground state energy 
per rung Aq, its O max is n by definition. 

For illustration, we consider the DES for the uniform 
spin ladder in second order, see Tab. [TTJ If we only target 
the ground state energy ho up to order 2, the maximum 



order of monomial A^ is given by 

°ma X (^4) = O max (,4o) - O m 

=*■ O raax (A 4 ) = 1. 



,(A 4 ) = 2-1 (33a) 
(33b) 



where we deal with equalities because there is only one 
contribution for dgho(l) in the DES and the O max (^4 ) is 
known. In this case the maximum order O max of A4 is 
lower than the targeted order 2. 

The O max of all coefficients can be calculated on the 
basis of the entire DES and of the minimum orders. Note 
that Eq. (32 1 defines O max implicitly, i.e., one has to 



find the correct self-consistent solution. This is done by 
starting from 



Omax (Ai 



if Ai is targeted 
otherwise. 



(34) 



A monomial Ai is targeted if we want to compute its co- 
efficient hi in the given order n. Starting from the initial 
choice (34) the Eq. (32) is iterated: O max is increased 
if necessary till convergence is reached. Convergence is 
guaranteed because we consider a finite set of {Ai} by 
construction and the O max (Ai) are bounded from above 
by n. Hence even in the worst case, there can be only a 
finite number of increments. For illustration, the maxi- 
mum orders for the uniform spin ladder in second order 
are given in Tab. [I] targeting dispersion or ground state 
energy. 

Once the maximum orders are known we can reduce 
the DES because some coefficients have a maximum order 
lower than their minimum order 



O max (^4i) < O m i n (Ai) 



(35) 



Thus they do not matter for the relevant quantities up 
to order n and can be discarded completely. Moreover, 
all contributions to the DES which use these terms can 
be neglected. In addition, all contribution can discarded 
for which 



O max (Ai) < O min {A 3 ) + O min {A k ) 



(36) 



holds. 

These considerations allow us to reduce the DES sig- 
nificantly. In Tab. [H] the reduction of the DES for the 
uniform spin ladder in second order is marked for the 
ground state energy (light gray) and for the dispersion 
(dark gray), respectively. 

We stress that one has to know the entire DES to apply 
the O max concept as described above. 



4- Simplification rules 

The reduction of the DES discards a large number of 
monomials and of the contributions Dij kl see for instance 
Fig. [3^,, which is essential for an efficient evaluation. But 
it would be even more advantageous if one avoided the 
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calculation of the omitted terms before they arc tediously 
computed. The minimum orders O min are known at each 
step of the iterative setup of the DES so that they can be 
used on the fly. But due to their implicit definition the 
maximum orders O max are not known during the set-up 
of the DES. 

Fortunately, estimates help. An upper bound for the 
maximum order is enough to accelerate the algorithm set- 
ting up the relevant part of the DES. Concomitantly, the 
memory consumption is reduced significantly^. Hence- 
forth, we call such estimates 'simplification rules'. Their 
concrete form depends on the structure of the perturbed 
and the unperturbed Hamiltonian, for instance the block- 
diagonality of the latter. We emphasize that the simpli- 
fication rules constitute the part of the epCUT method 
which depends on the model. 

In the following, we aim at a quantitative description 
up to order n of the effective Hamiltonian's blocks per- 
taining to at most q quasi-particles. For instance, q = 
provides the correct perturbative expansion of the ground 
state energy and q = 1 allows us to calculate the disper- 
sion relation up to order n. 

A monomial creating c triplons and annihilating a 
triplons is targeted if both c < q and a < q hold. Its max- 
imum order is the targeted order O max = n. If it is not 
targeted it can influence the targeted terms by affecting 
terms consisting of fewer creation and annihilation oper- 
ators via the DES. For the Heisenberg ladder, the unper- 
turbed Hamiltonian (Eq. (16)) is block-diagonal. Hence 
no commutation of generator terms with H changes the 
number of created and annihilated triplons. The leading 
order of the generator is 1, i.e., rj — O (x). 

In the commutation of a monomial with a generator 
term some of the local creation and annihilation opera- 
tors may cancel due to normal-ordering. In order to yield 
a term affecting the first subspaces with q quasi-particles 



c = max(c — q, 0) 
local creation operators and 

a' = max(a — q, 0) 



(37a) 



(37b) 



local annihilation operators have to cancel. 

First, we consider commutations with lowest-order 
generator terms stemming from the initial Hamiltonian. 
In the spin ladder, these terms have order 1 and create 
or annihilate AQP = 2 quasi-particles on adjacent rungs. 
Because each commutation with r/ 1 ) increases the order 
of the affected coefficients by one, the maximum order is 
bounded by 



~c'~ 




'a'' 


2 




~2 



(38) 



where the tilde on the left side means that one is deal- 
ing with an upper bound and \y\ stands for the smallest 
integer that is still larger or equal to y. If in the cal- 
culation of deH 1 -™' the estimate O max of a monomial is 



lower than m, this contribution is irrelevant and can be 
omitted. This reduces the size of both the DES and of 
the Hamiltonian to be tracked. Moreover, discarding ir- 
relevant monomials avoids the calculation of unnecessary 
commutators in the following iterations of the algorithm. 

Clearly, the number of created and annihilated quasi- 
particles can be reduced by a number AQP larger than 
2 by means of commutations with generator terms in- 
volving more quasi-particles which may have developed 
during the flow from the basic terms. But the gener- 
ator terms involving more quasi-particles have a higher 
minimum order O min so that a single commutation with 
them affects coefficients only in a higher order m + O min . 
In fact, for the used generator schemes, the ratio be- 
tween AQP and O min for new terms developed during the 
flow can not exceed the corresponding ratio for generator 
terms present in the initial Hamiltonian. Therefore, it is 
sufficient to consider only commutations with the initial 
terms in our simplification rules. 

The above generic simplification rule can be easily 
adapted to other models as long as the unperturbed 
Hamiltonian Hq is block-diagonal. Otherwise, Hq will 
lead to generator terms of order zero which means that 
terms with high quasi-particle number can influence the 
coefficients of terms with low quasi-particle number in 
the same order. This is why it is desirable to set up the 
perturbation in such a way that Hq is block-diagonal in 
the number of quasi-particles. 

Applying the simplification rule reduces the number of 
representatives considerably, see Fig. [3j leading to a sig- 
nificant improvement of runtime and memory consump- 
tion. This basic simplification rule can be improved fur- 
ther by taking more model-specific information into ac- 
count. A possibility to exploit the real-space structure 
of the monomials to lower the upper bound O max is de- 
scribed in Appendix [A} 

The computationally most costly part in the calcula- 
tion of the DES is the evaluation of commutators. Be- 
cause the simplification rules sketched above can only be 
applied after the commutation, we refer to them as a- 
posteriori rules. For the sake of efficiency, it is highly 
desirable to extend them to a-priori rules estimating 
whether a commutator has to be evaluated at all prior to 
its computation. We describe such a-priori simplification 
rules in Appendices [B] and [C| 

Because these rules are necessarily less strict than there 
a-posteriori analogues, one should use the combination 
of both kinds in practice. The additional use of a-priori 
rules does not reduce the number of representatives or 
the memory consumption. But it boosts the speed of 
the calculation significantly because the vast majority of 
commutators can be discarded, see Fig. [3Jd, and the a- 
priori rules help to avoid the laborious computation of 
these unnecessary commutators. 
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Figure 4. (Color online) Ground state energy per rung E of 
the symmetric spin ladder vs. relative leg coupling x in order 
17 using various evaluations. The direct evaluation (dark gray 
(blue)) renders a much more stable and reliable extrapolation 
of the plain perturbative series (black) than the various Pade 
extrapolations (light gray (orange)). 




Figure 5. (Color online) Ground state energy per rung Eo 
of the symmetric spin ladder vs. relative leg coupling x for 
different orders using the direct evaluation. 



F. Directly Evaluated epCUT 

In addition to the perturbative evaluation, the reduced 
DES computed by epCUT in a given order n can be 
evaluated non-pertubatively. After the reduction step 
described in Sec. |IIIE3| the DES consists exclusively of 
contributions which are relevant to the targeted quanti- 



ties in the desired order n. This reduced DES in Eq. (21a 



can be numerically integrated for any given value of x to 
obtain the coefficients of the Hamiltonian hi(£) directly 
without passing by an expansion in x. In such a cal- 
culation all coefficients influence one another to infinite 
order. The numerical solution depends on the expansion 
parameter in an intricate manner and can no longer be 



understood as finite partial sum of an infinite series. In 
this sense, the perturbative reduced DES in order n is ex- 
trapolated by the direct evaluation in a non-perturbative 
way. To stress the difference to perturbation series com- 
puted by epCUT, we call this technique directly evaluated 
epCUT (deepCUT). 

We emphasize that the reduction of the DES before 
the numeric integration is essential. It enhances the per- 
formance of the integration because the reduced DES is 
much smaller. But the crucial observation is that the re- 
duction renders the integration much more robust. Nu- 
merical integrations of the full DES diverge for high or- 
ders and high values of x. We conclude that the reduced 
DES represents the relevant physical processes in a more 
consistent way. The integration of the full DES generates 
spurious higher order contributions which overestimate 
certain effects. In an exact solution the spurious higher 
order contributions would be compensated by other pro- 
cesses which are captured only in a higher order calcula- 
tion. 

Analogous observations are known from diagrammatic 
perturbation theory where the inclusion of subsets of di- 
agrams in infinite order does not guarantee improved re- 
sults. Improved results can only be expected from sys- 
tematically controlled calculations. The inclusion of infi- 
nite orders is indicated if this achieves conserving self- 
consistent approximations. For instance, the shift of 
poles in a propagator is not captured by any finite pertu- 
bation series in the propagator, but it follows easily from 
a perturbation of the self-energy^. 

Figure [4] compares the ground state energy per rung 
i?o of the symmetric spin ladder as function of the rela- 
tive leg coupling x obtained from the plain perturbative 
series in order 17, from various Pade extrapolations, and 
from the direct evaluation. The plain series shoots up 
at about x ~ 0.7 while the Pade extrapolations start to 
scatter strongly beyond x ~ 1. The direct evaluation lies 
between the two stiffest Pade extrapolations and remains 
stable up to even very large values of the expansion pa- 
rameter ik3. Comparing various orders, see Fig. [5] the 
results of deepCUT converge rapidly and display only 
minor corrections for large values of x indicating a high 
reliability. In Sec. |IV| a comparison of the deepCUT re- 
sult with those of other methods will be presented. 

This deepCUT be ars similarities to the sCUT 
approacH^ 39 ' 4 " ' 52 ' 58 '. In sCUT, a set of basis opera- 
tors is selected by a truncation scheme and for this set 
the full DES is computed. It comprises all commuta- 
tion relations between the selected basis operators. In 
deepCUT, the order of the expansion parameter takes 
over the role of the truncation scheme. But we stress 
that deepCUT is not self-similar: In sCUT all commu- 
tators between the selected monomials are considered. 
In epCUT and thus in deepCUT only the commutators 
between specific subblocks based on the minimum or- 
ders O min are considered, see Sec. IIIC Moreover, tar- 
geting certain subspaces with q quasi-particles and the 
concomitant reduction of the DES does not only discard 
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irrelevant monomials. Also contributions linking relevant 
monomials are cancelled if their effect is of too high order. 
Therefore, the 'truncation' taking place in (de)epCUT, 
controlled by the expansion parameter, is a truncation of 
the DES rather than a truncation of operators as it is 
done in the sCUT approach. 

One practical advantage of the deepCUT over the 
sCUT is that only one parameter, the maximum order 
of the expansion parameter, needs to be fixed in order 
to define the approximation. In the sCUT, generically 
many parameters define the truncation scheme^EZESMMl 
which leaves some ambiguity about how to systematically 
improve the approximation. 

Another comparison of approaches is in order. Re- 
cently, Yang and Schmidt proposed a CUT approach 
based on graph theory (gCUTpSI. To compute a certain 
quantity such as the ground state energy the irreducible 
contributions of subgraphs, i.e., of linked clusters, of the 
lattice are summed. The size of the largest subgraph 
considered determines the approximation. The larger it 
is the better the system is described because physical pro- 
cesses with a larger range are kept. Thus the fundamental 
idea of the approach is similar to the one of deepCUT: 
Truncation in the range of processes, but local processes 
are kept to infinite order. 

The main difference is that the actual CUT is done 
on clusters. So gCUT has advantages and disadvantages. 
An advantage is that it is sufficient to deal with finite 
dimensional Hilbert spaces and the transformations can 
be performed on matrices. A disadvantage is that mo- 
mentum conservation cannot be exploited on the level of 
the clusters because they are finite which restricts the 
choice of applicable generators^. The deepCUT is based 
on second quantization^! and can take advantage of all 
symmetries of the problem under study. A detailed com- 
parison of both approaches is left to future studies. 



G. Transformation of Observables 



In order to calculate spectral densities for instance, 
the coefficients of the corresponding observable must be 
known with respect to the same basis as the effective 
Hamiltonian. Thus the observables must be transformed 
as well. This can be realized by integration of the flow 
equation for observables 



0(0) 0(1) 0(2) 0(3) Q(4) 



dtO(£) = [ri(l),0(£)} (39) 



introduced by Kehrein and Mielke^J^I. 

In analogy to the transformation of the Hamiltonian 
discussed in Sec. |III A| we introduce an operator basis 
Bi for the observable shifting the dependence on £ from 




Figure 6. Sketch of the epCUT algorithm to calculate the 
DES for the iterative calculation of de0^ 4 K Due to the com- 
mutators [n {1) ,O i3) ], . . . , [j? (4) , O (0) ], new terms with O min = 4 
emerge. In contrast to the algorithm for the Hamiltonian, see 
Fig. |2| no self-consistent calculation is needed for [r]^\ O^ ']. 
Self-consistency is required only for \rf ^ , O' 4 - 1 ] if is finite. 



the operators to their coefficients 
0(£) = Y.°^) B * 



EE/. 

i m— 



(m),obs//A m 



x m B l , 



(40a) 



(40b) 



where the second equation stands for the perturbative 
expansion of these coefficients. Hence the flow equation 
for observables ( 39 ) leads to a DES for their coefficients 



dM£) 



)o k {£). 



(41a) 



The contributions are obtained by calculating the 

commutators between the monomials of the generator 
and the monomials of the observable followed by a com- 
parison of the coefficients 



(41b) 



The differential equations ( 41a ) imply a hierarchical 
DES for the perturbative series ( |40b ) for the coefficients 

c^ m) ' obs m = E E D °3ktf\e)ft hohs (t)- (42) 

jk p+q=m 

The algorithm for the calculation of the DES in Sec. 
|III C| can be easily adapted for the transformation of ob- 
servables. Each order of the differential dtO^™"' is cal- 
culated recursively, cf. Fig. [6] Since the generator r\ is 
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defined solely by the Hamiltonian, it is not influenced by 
the outcome of the transformation of observables. For 
this reason the evaluation of [V m \ O*- -*] does not need 
to be carried out self-consistently. After the calculation of 
the commutators [77W, 0( m_1 )] . . . [r?^" 1 ), O^] , only 
the block [r)(°\ C^" 1 )] has to be treated self-consistently. 
But recall that r/°) only occurs if the unperturbed Hamil- 
tonian Hq is not (block-)diagonal. Because both differen- 
tial equations Q and (39) are coupled by the generator, 
their integrations have to be done simultaneously. 

For the transformation of the Hamiltonian we exten- 
sively discussed that only certain contributions really 
matter. We introduced the concept of a maximum or- 
der in which the coefficient of a physical process needs to 
be known in order to influence the targeted quantities. 
This concept allowed us to reach significantly higher or- 
ders. Thus we want to extend the concept of a maximum 
order also to the transformation of observables. It turns 
out that this extension is more subtle than one may think. 

Before, in the flow of the Hamiltonian, the maximum 
order of a generator coefficient O r ^ ax (Ai) is the maximum 
order of the same monomial 0^ ax {Ai) in the Hamilto- 
nian. Now we also target certain blocks of the observable 
and they are influenced by the monomials in the genera- 
tor. This leads to maximum orders for both the observ- 
able term 0^ ax {Bi) and the generator terms 0^; ax (Ai). 
The latter does not need to coincide with the maximum 
order 0"^ ax {Ai) resulting from the consideration of the 
Hamilonian flow alone. Thus one has to find a unique 
and unambiguous way to fix 0^ ax (Ai). We discuss three 
alternatives: 

(A) The maximum order of the generator terms is cho- 
sen in such a way that the targeted quantities in both the 
Hamiltonian and the observable(s) can be computed up 
to the targeted order^ n 

OlM) = maxCO^A-), 0«&(Ai)). (43) 

Then the iterative calculation of the O max must be real- 
ized within a single self-consistent loop. The perturbative 
evaluation yields a perturbative series for the coefficients 
of the observables under the transformation with the full 
generator up to order n. It may happen that in this way 
some generator terms are assigned a higher 0^ ax > OHi ax 
than in the transformation of the Hamiltonian alone so 
that the DES of the Hamiltonian comprises additional 
contributions. By construction, this does not affect the 
perturbative evaluation of the epCUT. But it will affect 
its direct evaluation (deepCUT) although it should be 
absolutely minor in a parameter regime of good conver- 
gence of the flow. 

(B) Alternatively, the determination of 0^ ax (Bi) and 
0^i ax (Ai) can be realized after and strictly separated 
from the calculation of 0^(Ai) and 0^ ax {Ai). Mono- 
mials which are discarded due to the reduction of the 
Hamiltonian will not be considered for the DES of the 
observables even though this may affect the targeted co- 
efficients of the observable. Hence the transformation of 
the observables in perturbative evaluation is not realized 



with respect to the complete generator. We stress that 
this does not violate the unitarity of the transformation 
up to the order of calculation since the generator is still 
antihermitian and it is essentially the same as for the 
transformation of the Hamiltonian. No signficant devia- 
tions are expected in the regime of good convergence of 
the flow. Note also that any generator whose coefficients 
differ only by orders larger than 0^ ax (Ai) leads to the 
same perturbative series for the relevant quantities in the 
Hamiltonian. 

(C) A third alternative consists in taking over the 
0^[ ax (Ai) for the reduction of the DES for the observ- 
ables. Then only the values 0° ax (Bi) are computed self- 
consistently. 

For deepCUT, alternatives (B) and (C) ensure that 
the DES for the Hamiltonian is independent of the con- 
sidered observables. Generally, we expect that the pre- 
cision in the derivation of effective Hamiltonians is more 
important than the precision of matrix elements. Also 
in experiment, energies are generically known to much 
higher accuracy than matrix elements. 

In order to keep the effective Hamiltonian in direct 
evaluation independent of the observables, we decide to 
use alternative (B) for deepCUT. For the perturbative 
evaluation, however, we favor alternative (A) because it 
makes the rigorous determination of the perturbation se- 
ries of matrix elements possible. 

The computational performance can again be increased 
decisively by applying simplification rules. They can be 
used directly for observables if both the Hamiltonian and 
the observables meet their requirements. This can often 
be achieved by appropriate definitions. For instance, the 
observable 



25 ' — t q + t z ^ Q + it t x l 



Kfityfi (44) 



is needed for the calculation of the dynamic structure fac- 
tor relevant for inelastic neutron scattering. But this ob- 
servable includes non-block-diagonal terms in order zero. 
To circumvent this problem, we consider the observable 
x ■ Sq' z instead. In this way, the non-block-diagonal 
monomials in the observable are shifted to order 1 so 
that they behave like the non-block-diagonal pertuba- 
tion in the Hamiltonian. One loses an order of accuracy 
for a given fixed order n of the calculation. But all the 
simplification rules relying on block-diagonality in order 
zero can be used as before. 



IV. RESULTS FOR UNIFORM SPIN LADDER 

A. Ground State Energy 

The ground state energy per rung is calculated up to 
and including order 17. The results of the direct (black 
solid line) and of the perturbative evaluation (dashed 



black line) are displayed in Fig. [7 



the perturbative series, see Tab. IV 



The coefficients of 
agree perfectly with 
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Figure 7. (Color online) Ground state energy _Eo per rung 
vs. relative leg coupling x resulting from various methods. 
The epCUT results (order 17; direct: black, solid; perturba- 
tive: black, dashed; Pade[ll,6]: black, dotted) agree with the 
sCUT results (d=(12, 10, 10,6,6, 5, 5,4,4), light gray (orange), 
solid) and the DMRG results (dark gray (blue)). The ener- 
gies from the direct evaluation, the sCUT and the DMRG 
lie on top of each other, see also upper inset. The deviation 
\AEo\ = \Eq, direct — -Eo,dmrg| is shown in the lower inset. 

the fractions from pCUT (up to order 14^1 and with 
the decimal numbers (up to order 23) given by Zheng et 
alP^l. This agreement shows that the epCUT works for a 
system with equidistant spectrum in Hq. 

The plain series is trustworthy up to x ~ 0.7. For 
larger x extrapolations are needed. The dotted black 
line shows the "best", i.e., stiffest, Pade-extrapolant of 
order [11,6] for this series. Other Pade-extrapolants are 
shown in Fig. [4] The result of the deepCUT are de- 
picted as solid black line in Fig. [7] It fits perfectly to 
the perturbative result for weak leg couplings. For larger 
values of a;, it serves as an excellent extrapolation of the 
perturbative results. In order to support that the di- 
rectly evaluated results are quantitatively reliable even 
for larger x the deepCUT data is compared with results 
from sCUI^Hl and from DMRCPHSfi] Tll(3 sCTJT re _ 

suits, represented by the solid light gray (orange) line, 
are calculated with the ground state generator and the 
truncation scheme d=(12,10,10,6,6,5,5,4,4). The trunca- 
tion reads d = {d2,d%, ^4, • ■ ■ ), where di denotes the real 
space extension of a monomial with i interacting quasi- 
particles. A monomial with i interacting quasi-particles 
is truncated if it exceeds the extension di. For instance, 
the monomial t* a r t a r+4 has an extension d2 = 4. 

The DMRG data, represented by a solid dark gray 
(blue) line, results from a finite-size scaling. The ground 
state energies for ladders with L = 40, 60, . . . , 160 rungs 
and m = 500 states are extrapolated with the ansatz 

e -L/L 

Eo(L) =SoM + co-^7- (45) 
to estimate the ground state energy for an infinity ladder 




0.0 0.2 0.4 0.6 0.8 1.0 

k[jt] 



Figure 8. (Color online) Dispersion u>(k) for various values 
of a; £ {0.5,0.8,1,1.5,2} and various evaluation techniques 
based on order 15. At k = 7r /2, the lowest curve is x — 0.5 
and the highest curve is x — 2. For x = 0.5, the plain series in 
x (light gray (orange), dashed) and for x = 0.8 and x = 1 the 
plain series in the parameter p(x) from Eq. ( |47[ ) (dark gray 
(blue), dashed) are shown. The results of the direct evaluation 
(black, solid) agree well with the perturbative results for small 
x and they are still robust for larger x. 



Eq(oo) for each value of x. 

The results of these methods agree perfectly. The devi- 
ations between the DMRG results and the results of the 
direct evaluation are shown in the lower inset of Fig. [7] 
They increase with rising x, but they remain still small. 
For x — 1.5, the deviation is less than 10~ 3 Jj_ and for 
x = 3 it is still less than 10~ 2 Ji. 



B. Dispersion 

The one-triplon dispersion is calculated up to order 15. 
The dispersion is obtained by a Fourier transform of the 
one-triplon sector of H c r . The dispersion reads 

n 

uj(k,x) = t + ^2t d cos(dfc) , (46) 

d=l 

where td is a hopping element over the distance d. Figure 
[8] shows the dispersion for various values of x. The plain 
series of the perturbative evaluation is depicted as dashed 
light gray (orange) line. The coefficients of this series 
agree quantitatively with other series expansion results^ 
up to order 8. 

The coefficients of the perturbative series of the spin 
gap match those of other series expansion methods^! up 
to order 13 . The plain series is reliable up to x « 0.6. 
For larger values of x extrapolations are needed. For the 
dispersion we used the extrapolation scheme based on a 
re-expansion of the original series in terms of a suitable 
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Figure 9. (Color online) Gap A(x) vs. relative leg coupling 
x for various orders (6, 10, 14, 15) using the direct evaluation 
(black, solid), perturbative evaluation (plain series; order 15; 
black, dashed) and a 1 /l 2 finite-size scaling DMRG result 
(dark gray (blue), solid). The lowest curve in direct evalua- 
tion at x — 3 stems from order 6; the highest curve from order 
15. The 15 th order curve in direct evaluation agrees very well 
ith the DMRG results. The deviations to the DMRG results 
are depicted in the inset. 



internal parametei^p(a;) 



p(x) = 1 - 



A(x) 
(l + x)J± 



(47) 



where A(a;) denotes the gap. In order to use the above 
mapping i-ipa reliable extrapolation of the gap A(x) 
is needed. This was achieved by dlog-Pade extrapolations 
which work very robustly for the gap. The dashed dark 
gray (blue) lines in Fig.|8]represent the plain series in this 
internal parameter without any further extrapolation. 

The solid black lines are the results of the direct eval- 
uation. For small x the perturbative and the direct eval- 
uation agree very well. The direct results are again very 
robust for larger x as well. 

To corroborate that the deepCUT results are reliable 
even for relatively large values of x its spin gap is com- 
pared to the one obtained in DMRdsMM - m pig. [g| The 
solid black lines represent the deepCUT results. The 
solid dark gray (blue) line depicts the results of a 1 /l 2 
finite-size scaling of the DMRG results based on the 
ansatz 



A (z) — A(oo) + a\z + a,2Z 3 / 2 + a%z 



(48) 



with z :— 1 /l 2 . 

The deepCUT data shows that for larger x a higher 
order is needed to compute the gap accurately. This is 
understood on the basis of the correlation length of the 
system. A larger value of x implies a larger correlation 
length £ which is given by v /a where v is the spin wave 
velocity if no gap were present. By construction of the 
(de)epCUT, the order n defines the range of processes 



which are still captured. Hence one can expect a reliable 
result as long as 



> 



> 



>/a(x), 



(49) 



where the lattice constant is set to unity. The velocity v 
can be estimated by fitting v sin(fc) to the maximum of 
the dispersion. We find indeed that ( 49 ) is satisfied up to 
x m 3 for n — 15. The deepCUT curve for n = 15 agrees 
well with the DMRG-results. The deviations shown in 
the inset are rather small. For x = 2, it is below 1CP 2 J±. 
Furthermore, the dispersions shown in Fig. [H] agree with 
exact diagonalization results^. 

We conclude that the reliable results for the uniform 
spin ladders beyond x — 1 illustrate the efficiency of the 
deepCUT approach for a model with an equidistant spec- 
trum. 



C. Spectral weights 

Here we use the transformation of self-adjoint observ- 
ables 



O (cf. Sec. IIIG) to address the issue of spec- 
tral weights. We denote the subspace spanned by the 
states with q quasi-particles by QP q . As in previous 
wor ]j26i2 7 | 30 | 56 | we S pij^ ^he total spectral weight at zero 
temperature into its contributions from the different sub- 
spaces QP q 



£ k*i°i°>i 2 

|i)eQP„ 



= <o|(M|ci) 



(50a) 
(50b) 



where 0| stands for the sum over all terms of the trans- 
formed observable consisting of q creation operators and 
p annihilation operators in normal-ordering. The state 
|0) denotes the vacuum state of the effective model, i.e., 
the ground state of the Hamiltonian. 

If the subspaces QP g have been separated by the CUT, 
i.e., the effective Hamiltonian does no longer mix them, 
the spectral weights d efined by (50) coincide with the 
ones defined previousl y ^ 26 * 27 !. The spectral weights corre- 
spond to the integral over momentum and and frequency 
of the corresponding dynamic structure factor S q (k.uj) 
where the subscript q denotes the contribution of the 
subspace QP g . Thus, separate sum rules exist for each 
QPq. Such a split- up is only possible because the dynam- 
ics does not mix the subspaces according to the above 
assumption. We recall that the dynamic structure fac- 
tors encode the response of various inelastic scattering 
experiments. 

If the subspaces QP are not or not all separated the 



equal-time definition (50) is still well-defined. But /„ 



can no longer be interpreted as the sum rule of S q (k, u>) 
because the subspaces mix in course of the dynamics in- 
duced by the Hamiltonian. Nevertheless, the values I q 
provide a plausible measure of the importance of the sub- 
spaces of different number of excitations. A large spec- 
tral weight for low numbers of quasi-particles indicates 
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Figure 10. 

(Color online) Spectral weights for different observables 



defined in Eq. (51 1 and numbers of triplons vs. relative leg 
coupling x. Panel (a) shows the S = 1 observable O , panel 
(b) and (c) show the parallel and perpendicular S — 
observables O 11 and O m . The calculations were carried out 
to order 8 for the modified observables x ■ O 1 , x ■ O , and 
x ■ O™ using deepCUT with the generator scheme 7/2- 
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that results of scattering experiments can be understood 
from the spectral densities involving only low numbers of 
quasi-particles. 

In the present work, we concentrate on the spectral 
weights for the observables 

(51a) 
(51b) 
(51c) 

to illustrate the transformation of observables. The 
observable O 1 induces a local spin S — 1 excitation 
which ca n be stud ied experimentally by inelastic neutron 
scatterin g 28 1 29 1 44 1 . The observables O 11 and O m induce 
5 = excitations which can be studied by optica l prob es, 
e.g., Raman scattering^ or infrared absorptio n 1 24 * 43 !, in 
polarizations parallel and perpendicular to the ladder, 
respectively. Because triplons are 5 = 1 states, both 
observables O 11 and O lu induce no contributions in the 
one-triplon channel: 1\ — 0. The calculation of the corre- 
sponding spectral densities S q (k, uj) is left to future stud- 
ies. 

Since the description in terms of triplons on rungs is 
obviously best for low values of x we expect that more 
and more triplons need to be addressed upon increasing 
x. To assess the relative importance of different triplon 



channels we introduce the relative weights I q . rc \ = I i/itot- 
They can be calculated using the sum rule 

oo 

hot :=^/ 9 = (0|OO|0)-(0|O|0) 2 (52) 

for the total spectral weight /tot- For the observables 
defined in (51 ), the total weights are given by 



1 tot — 
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with the variables 
Y :=2(0|S 

Z:= (0|S 



sf |0> 



dx 



So |0) = E 



dE 
dx 



(53a) 
(53b) 
(53c) 

(54a) 
(54b) 



where we use the ground state energy per rung Eq. 

We focus on the spectral weights in the first four- 
triplon (four-quasi-particle) channels I\, 72,^3 , and I4 up 
to large values of the relative leg coupling x = 3 us- 
ing deepCUT. In this region, a complete decoupling of 
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all subspaces using ?7p C or 774 is no longer possible be- 
cause divergences occur in the numerical integration of 
the flow. This problem is well-known from sCUT; it 
stems from th e overla p of continua of different quasi- 
particle numbe d * 51 * 52 !. In this situation, the sorting of 
quasi-particle spaces ascending by energy is no longer 
possible. In the perturbative evaluation of epCUT, no di- 
vergencies appear if the energies in Hq are separated and 
indeed ordered according to ascending number of quasi- 
particles because the hierarchy in Eq. ( 26 1 precludes any 



feedback of high order coeffients to low order coefficients. 
In the alternating spin ladder the epCUT based on the 
quasiparticle number must be modified for y > 3. But 
we stress that this reflects a more sophisticated physics 
which must be considered in the choice of the generator. 
It does not represent a conceptual problem of epCUT. 

To avoid convergence problems due to overlapping 
continua in deepCUT, we aim only to decouple sub- 
spaces with at most 2 quasi-particles using r]2 while keep- 
ing monomials linking subspaces with higher number of 
quasi-particles, for instance QP 3 O QPs- As a con- 
sequence, the observables are transformed to a quasi- 
particle basis where 3 and 4 quasi-particle states still 
couple to other subspaces. 

The spectral weights for the S = 1 observable O 1 
depict ed in p anel (a) of Fig. 10 agree well with pCUT 
result^* 27 * 53 * for small values of x. Note that only the 
one- and the two-quasi-particle channel can be com- 
pared quantitatively because the pCUT separated also 
the three- and four-quasi-particle subspaces, but the 
present calculation does not. 

For O l , most weight is concentrated in the first two 
quasi-particle channels. Even at x = 3, the one-triplon 
channel still contains 57.9% of the total weight. The rel- 
ative weight of the two- and three-triplon channel rises 
up to 35.2% and 11.3%, respectively. The four-triplon 
weight remains negligible. The sum rule is slightly vi- 
olated because the accumulated relative weights exceed 
100%. This inaccuracy is related to the finite order of 
calculation. The degree of the violation of the sum rule 
can be used as a measure for the reliability of the results. 
Even at x = 3, the excess weight is only 5.5%. 

Panel (b) of Fig. 10 shows the spectral weights I q for 
the S = observableD 11 . Since triplons have spin S = 1, 



there cannot be any weight in the one-triplon channel. 
Instead, most weight is concentrated in the t wo-triplon 
channel which agrees well with pCUT results^ 21 - 26 1 27 1 53 1 . 
Compared to O l , the three-triplon channel is much more 
pronounced displaying a relative weight of 44.4% at x = 
3. At this value, the sum rule is fullfillcd within 6.7%. 

The observable O lu is symmetric with respect to the 
ladder's centerline. This impli es tha t this observable does 
not change the parity of a stat d 26 * 27 -l A single triplon is an 
odd excitation with respect to the ground state. Hence 
O m can create or annihilate triplons only in pairs. There 
is no weight in odd channels in Fig. [lO^c). As a conse- 
quence, the spectral weight is distributed over the two- 
and four-triplon channels only. Our results do not indi- 
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Figure 11. (Color online) Ground state energy per rung 
Eo vs. relative leg coupling x for various values of y £ 
{1,1.1,1.2,1.3,1.4,1.5} and various evaluation techniques 
based on the DES in order 16. The highest curve at x = 
is y = 1 and the lowest curve is y = 1.5. Again the direct 
evaluation (dark gray (blue)) yieldsw a much more stable and 
reliable extrapolation of the plain series (black) than the var- 
ious Pade-extrapolants (orange). 



cate a sizeable contribution from six and more triplons, 
but this has not been studied quantitatively. At x = 3, 
the sum rule is violated by 7.5%. Indeed, this violation 
sets in at about x — 0.6 when the four-triplon weight 
becomes significant. Thus we presume that the latter is 
a bit overestimated, but we could not identify the mech- 
anism for this effect. The perturbative results for the 
weights fulfill the sum rule to the required order. 



V. RESULTS FOR ALTERNATING SPIN 
LADDER 

A. Ground State Energy 

For Jj_ J^, the ground state energy is calculated 
up to order 16. Due to the doubled unit cell, only a 
slightly lower order can be reached than for the uniform 
spin ladder. Roughly, we need double the number of 
coefficients for the alternating spin ladder. The ground 
state energy per rung is given by Eq — h o/2. 

The perturbative results from epCUT are shown in Fig. 
[TT] As expected the ground state energy decreases upon 
rising y. The black lines represent the results of the plain 
series for various y. The coefficients are given in Tab. |IV| 
The light gray (orange) lines correspond to various Pade- 
extrapolants. The plain series is reliable up to x ~ 0.75 
for y = 1 and up to x ps 0.85 for y = 1.5. So the x 
up to which the series is reliable depend on the value 
of y. Since a larger value of y supports the dominance 
of the unperturbed Hamiltonian Hq it is clear that an 
increasing y supports the validity of the perturbation. 
The dark gray (blue) lines represent the results of the 
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Figure 12. (Color online) Ground state energy Eq per rung vs. 
relative leg coupling x for y — 1.2 for various methods. De- 
picted are the results of the perturbative evaluation (order 16; 
black, solid) the direct evaluation (order 16; black, dashed), 
a high-level sCUT calculation (d=(12,10,10,6,6,5,5,4,4), light 
gray (orange), solid) and a DMRG calculation (dark gray 
(blue), solid). The direct evaluation agrees very well with 
the sCUT and the DMRG results. The deviation \AE \ = 
| -^O, direct — £?o,dmrg| between the results of the DMRG and 
of the direct evaluation is shown in the lower inset. 



deepCUT. These results again represent a very robust 
extrapolation of the perturbative results up to larger x. 

To show the efficiency of the epCUT the results for 
the ground state energy per rung are compared to the 
results of an sCUT calculation and a DMRCpHl!] calcu- 
lation. The sCUT was performed with the ground state 
generator and a d=(12, 10, 10,6,6,5,5,4,4) truncation. The 
DMRG results (L = 20, 40, ... , 100, m = 100) are ex- 
trapolated according to Eq. (45). Figure [12] compares 
the results of the various approaches. They agree very 
well with one other. The deviations between the results 



of the direct evaluation and the DMRG calculation are 
small, see upper inset, are small. For x — 1.5 the devia- 
tion is less than 10~ 3 Ji . 



Figure |l3[a) displays the dispersions for x = 0.5 and 
various values of y. The solid lines represent y = 1, the 
dashed ones y = 1.2 and the dotted ones y = 1.4. The 
dark gray (blue) lines stand for directly evaluated results 
and the black lines for the perturbatively evaluated ones. 
For x — 0.5 the plain series is used. Both results agree 
very well. For y = 1 we retrieve the uniform ladder 
and the two branches meet at k = 1/2. As expected the 
branches split at k = n '/ '2 once y > 1 holds due to the 
reduced translational symmetry. 

To show the efficiency of the epCUT, the disper- 
sion relations for y — 1.2 is compared to the disper- 
sion from an sCUT calculation in Fig. [II] The sCUT 
was performed with the generator rji and the trunca- 
tion ci=(12,10,10,6,6,5,5,4,4). The dispersions match per- 
fectly. The deviation between sCUT and the perturba- 
tive evaluation is less than 10 _3 J_l. The differences in 
the upper branch are larger than the ones in the lower 
branch. 

Furthermore, the gap for y = 1.2 is compared to the 
gap obtained by a DMRG calculatiorPHSl. The finite- 
size scaling is carried out again based on Eq. (48). In 
Fig. [15] the solid black line shows the perturbative result 
and the dashed black line the result of the deepCUT. The 
result of the DMRG calculation is depicted as solid dark 
gray (blue) line. The deviations between the deepCUT 
and the DMRG calculation are shown in the inset. Again 
the results agree very well, e.g., the deviation is less than 
10~ 2 J± even at x = 2 . 

The dispersions at x = 1 of the direct and of the per- 
turbative evaluation are plotted in Fig. 13 b). The per- 



turbative result are rendered using the plain series in an 
internal parameter^. The parameter, however, defined 
in Eq. (47) does not work because at y 7^ 1 it behaves 
like p(x) (x x 2 and not linearly in x anymore. Thus the 
series in x cannot be re-expressed in a series in p. So we 
modify the internal parameter 



(M e , 



M n , 



2|M H , 



(56) 



B. Dispersion 

The dispersion is calculated up to order 13 for the al- 
ternating ladder. An important step is the Fourier trans- 
form of the hopping in the one-triplon sector of H e g. But 
the doubled unit cell has to be taken into account which 
halves the Brillouin zone. In return, the dispersion ac- 
quires two branches reading 



w±(fc) 



M eo + M 00± J(M ee -M oo r +ML 



where stands for the Fourier transform of the hop- 
ping processes from a rung of parity i to a rung of parity 



where all matrix elements are taken at vanishing wave 
vector k — 0. We choose this parameter because it re- 
produces the previous definition (47) for y = 1 and for 
x — > 00. In addition, it fulfills p a cx x for x — > for all val- 
ues of y. Otherwise, the extrapolation can be performed 
as before^. The Fourier transformed matrix elements 
Mij(x) are obtained by robust dlog-Pade extrapolations. 
The results of deepCUT and of the series in this internal 
parameter agree very well. 

The epCUT results for the alternating ladder exem- 
plify the efficiency of this CUT for a system with a non- 
equidistant spectrum in Hq . Thereby, the range of appli- 
cability of perturbation by CUT s is crucially enhanced 
because the previous pCUT^^U is restricted to equidis- 
tant unperturbed spectra. 
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Figure 13. (Color online) Panel (a): Dispersion uj(k) for x = 0.5 and various values of y £ {1, 1.2, 1.4} (solid, dashed, dotted; 
order 13). The direct evaluation is depicted by dark gray (blue) lines and the perturbative one by black lines which are hardly 
visible because they are just below the other lines. For the perturbative results the plain series are used. Panel (b): Dispersion 
oj(k) for x — 1 and various values of y £ {1, 1.2, 1.4} (solid, dashed, dotted; order 13). The direct evaluation is depicted by dark 
gray (blue) lines and the perturbative one by black lines. For the perturbative results the plain series in the internal parameter 
p a (x) defined in Eq. ( |56[ ) without any further extrapolation is used. The upper branches differ slightly, but in general both 
evaluation techniques agree well. 




0.0 0.1 0.2 0.3 0.4 0.5 

k[%] 



Figure 14. (Color online) Dispersion w(fc) for x = 0.5 and 
y — 1.2. Displayed are the direct evaluation (order 13; 
black, dashed) and the perturbative one (order 13; black, 
solid), and an sCUT result (d=(12, 10,10,6,6,5,5,4,4), dark 
gray (blue), solid). All dispersions lie on the top of each other. 
The deviations |A u /jj_| = \ u &vtect/j ± — ^pxn/ (black) and 
|A w /j_,_| = | w »cuT/jj_ — "pmt/jjj (blue(dark gray)) are shown 
in the inset. The solid (dashed) lines depict the lower (higher) 
branch. 

VI. CONCLUSIONS 

In this article, we presented a methodical develop- 
ment and illustrated it for a well understood model. We 
extended the previously known perturbative continuous 
unitary transformation (pCUT) in two ways. 

First, we formulated the perturbative realization of the 
CUTs directly in second quantization. Thereby, the un- 
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Figure 15. (Color online) Gap A(x) vs. relative leg coupling x 
for y — 1.2 from the pertubative evaluation (order 13; black, 
solid), the direct evaluation (order 13; black, dashed) and a 
DMRG result (dark gray (blue), solid) extrapolated to the 
thermodynamic limit L = oo by finite-size scaling oc 1 /l 2 . 
The results agree well. The deviations between deepCUT 
and DMRG results are plotted in the inset. 

perturbed part is no longer restricted to an equidistant 
spectrum of energy eigen values. The direct expansion 
of all coefficients in the effective Hamiltonian is not effi- 
cient enough. But tracking the powers in the expansion 
parameter x of all the physical processes it is possible 
to identify the relevant ones for the low-energy effective 
model: Ground state energy, single quasi-particle disper- 
sion, and two-quasi-particle interactions. We could show 
that this leads to an efficient and competitive approach to 
obtain effective models. Their parameters are computed 
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as series in the expansion parameter. For distinction, we 
baptized the enhanced approach enhanced perturbative 
CUT (cpCUT). 

Second, we found that the system of differential flow 
equations, which has been reduced to provide the per- 
turbative series representation of the effective model, can 
also be directly evaluated. It appears that this directly 
evaluated perturbative CUT (deepCUT) yields a very ro- 
bust and reliable way to exploit the information in the 
perturbative differential flow equations. In some sense, 
one can think of it as a robust extrapolation though we 
stress that it is not an algorithm applied to a series. The 
deepCUT provides the parameters of the effective models 
for given initial Hamiltonian. Each set of initial parame- 
ters requires a numerical integration of the flow equations 
which is a moderate numerical task. The essential effort 
lies in deriving the system of differential flow equations 
which is the same as for the epCUT. 

The two abstract key results were illustrated by calcu- 
lations for antiferromagnetic S — 1/2 spin ladders with 
two legs. Two types of spin ladders were studied. The 
expansion parameter is the leg coupling relative to the 
(smallest) rung coupling. The uniform spin ladder with 
the same value of the rung coupling is the standard model 
which is very well studied. The alternating spin ladder 
with alternating rung couplings has not yet been studied 
to our knowledge. For the present purposes, it constitutes 
a model with a non-equidistant unperturbed spectrum if 
the perturbation is set up around the rung Hamiltonian. 

For the uniform spin ladder, the known series coeffi- 
cients could be retrieved by epCUT. The corresponding 
results for the alternating spin ladder have not been pub- 
lished elsewhere. They show that general unperturbed 
spectra can be treated by epCUT. 

The data obtained by deepCUT illustrate that this ap- 
proach yields surprisingly robust results. The uniform 
spin ladder could be treated up into the strong-leg limit 
with values of the relative leg coupling x = Jn/J± of up 
to x = 3. This is a par amet er regime which was not 
accessible by CUTs beforcPHl. 

The limit of the applicability of deepCUT can be un- 
derstood in terms of the correlation length. A deepCUT 
calculation in order n in a perturbation linking adjacent 
sites allows us to capture processes up to the range n ■ a 
where a is the lattice constant. Hence reliable results can 
be expected if the correlation length £ = v/A is lower 
than n ■ a. The deepCUT results for the alternating lad- 
der are also very robust, though a little less than for the 
uniform ladder. 

The idea behind the deepCUT approach bears some 
similarity to the graph-based CUT (gCUT) recently in- 
troduced by Yang and SchmidtP^. In the gCUT the ir- 
reducible contributions of subclusters of increasing num- 
ber M of sites of the lattice are summed. The number 
M plays a similar role as the number n in deepCUT. It 
restricts the maximum range of the tracked physical pro- 
cesses. These processes are treated in infinite order in 
the ratios of the coupling so that the evaluation is non- 



perturbative. Noticeable differences between deepCUT 
and gCUT are the following. 

The gCUT approach performs calculation on finite 
clusters and can be done on the level of states in 
the finite-dimensional Hilbert space, assuming that the 
Hilbert space of a single site is finite. Only the Hilbert 
space of the clusters with up to M sites needs to be con- 
sidered. Translational invariance is lost on the level of 
the single cluster which limits the choice of applicable 
generators. 

The deepCUT approach works on the level of monomi- 
als of creation and annihilation operators, i.e., in second 
quantization. Thus essentially all symmetries of the lat- 
tice problem under study can be preserved by construc- 
tion. A large variety of the generators can be realized. A 
detailed comparison of deepCUT and gCUT is an inter- 
esting issue left for future research. 

In a nutshell, we advocate two approaches to derive 
effective models in a systematically controlled way in this 
article. They have been illustrated for spin ladders and 
we expect that applications to many other models will 
soon be possible. 
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Appendix A: Extended A-Posteriori Simplification 
Rule 



The upper bound O max for the maximum order can 
be reduced by considering the real-space structure of 
the monomial. For clarity, we restrict ourselves to one- 
dimensional models. As for the basic simplification rule, 
we discuss the effect of commutations with first-order 
terms present in the initial Hamiltonian. This is sufficient 
because any more complicated monomials in the genera- 
tor have been induced by commutations of a number of 
first-order terms. Hence their gain in number of involved 
quasi-particles is paid for by a correspondingly higher or- 
der in x. Thus one may safely restrict the consideration 
to the basic building blocks present in the initial Hamil- 
tonian. 

For the spin ladder in terms of triplon operators (Eq. 



(9d)), a commutation with the generator rf 1 ' cancels at 
most two local creation or annihilation operators on ad- 
jacent sites. Therefore, sparse and extended monomials 
require more commutations in order to reduce their local 
operators compared to monomials with the same num- 
bers of operators which are more localized in real space. 
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Figure 16. Decomposition of the sites with creation operators (or the annihilation operators, respectively) of a monomial into 
linked subclusters ki and its covering with first-order generator terms. Each circle stands for a rung of the spin ladder. Filled 
circles represent rungs where the local action of the monomial differs from identity. At most two adjacent local operators can 
be cancelled by a single commutation with , this is represented by ellipses. 



At first, we study the ground state energy per rung, 
i.e., the coefficient of the identity operators summed over 
all rungs, in highest order. The clusters of the creation 
and of the annihilation operators are treated separately. 
Both are decomposed into linked subclusters of size kf 
and kf, see Fig. 



1G 



To cancel all local operators, each 
subcluster needs to be covered by [%1 first-order gener- 
ator terms. In conclusion, 



(Al) 



commutations with 7/ ' are needed for the clusters of 
creation or annihilation operators to be reduced to the 
coefficient of the identity operator. This argument leads 
to the extended upper bound for the maximum order 



On 



= n - K£ - K a . 



(A2) 



For a single linked cluster, this formula resembles the 



result obtained for the basic simplification rule ( 38 1 



The formula ( A2 ) can be generalized to 



O max = n - K c „ - K" 



(A3) 



for the subspace QP ? of states with q quasi-particles lead- 
ing to modified cluster sums K q . Let q be the number of 
the targeted subspace with the highest number of quasi- 
particles. This means that q ist the maximum number 
of local creation and annihilation operators allowed in a 
monomial targeted up to order n. Terms which affect 
more quasi-particles have to be reduced to affecting at 
most q quasi-particles by commutations with 77W until 
at most q local creation and annihilation operators are 
left. To obtain an upper bound O max , one has to choose 
q positions for local operators to be kept in the cluster in 
such a way that the other creation and/or annihilation 
operators can be cancelled by a minimum number of com- 
mutations. To this end, one also has to consider that the 
commutations with hopping terms stemming from 
may also shift creation and/or annihilation operators so 
that they form adjacent pairs which can be cancelled by 
pair creation or annihilation. But it turns out that this 
mechanism can reduce the cluster sum Kq at most by 
unity while the elimination of a pair of adjacent local op- 
erators always reduces the cluster sum by unity. Hence 
the latter process dominates and provides the correct up- 
per bound Omax- 



For the hopping in the symmetric ladder model, the 
above approach means to select sites at the edges of odd 
subclusters first. This saves one commutation for each 
local operator kept. Let a be the number of odd clusters. 
The cluster sum Kq is reduced in this way by 



d\ = min(a, q). 



(A4) 



If more local operators remain, i.e., a < q, the most 
efficient way to place them is in pairs on even subclusters. 
This reduces the cluster sum additionally by 



di 



(A5) 



where [y\ is the largest integer which is still smaller or 
equal to y. 

In conclusion, the cluster sums are reduced when one 
is aiming at higher quasi-particle subspaces according to 



K' q =K -d x - 



q-di 



q + di 



(A6a) 



To avoid unreasonable negative results, this expression 
has to be checked against zero to obtain the final result 



K q = max(K'O). 



(A6b) 



We remark that the extended simplification rule can be 
easily adapted to other models with monomials of first 
order in the generator create or annihilate an arbitrary 
number AQP of quasi-particles on adjacent sites. For 
further refinements of O max , one may the triplon polar- 
izations x, y, z as well. But the derivation and applica- 
tion of an appropriate polarization-sensitive simplifica- 
tion rule is beyond the scope of this article which aims 
primarily at the proof-of-principle demonstration of ep- 
CUT and deepCUT. 



Appendix B: Basic A-Priori Simplification Rule 



As stated in Sec. HIE 4 the performance of the epCUT 
algorithm can be enhanced significantly by avoiding the 
computation of unnecessary commutators. For this pur- 
pose, we consider the two normal-ordered products TD 
and DT in 



[T, D] = TD - DT 



(Bl) 
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separately. Here we discuss TD explicitly; DT is treated 
in the same way. For an analogue of the basic sim- 
plification rule (Sec. HIE 4), we estimate the minimum 
numbers of creation and annihilation operators ctd and 
o-td which can appear in the monomials of the normal- 
ordering of TD. We use the numbers ct , cd , ax , and ajj 
from each factor as input. At most 



s TD = mm(a T ,c D ) 



(B2) 



pairs of local operators can cancel in the process of 
normal-ordering. Hence it follows 



c T D > c T + c D - s T D 
a T D > ax + a D - std- 



(B3a) 
(B3b) 



Using these estimates in Eq. ( 38 ), one obtains an upper 
bound 



O 



max,TD 



max 



ct + c D - s T D 



aT 



2 



-q,0 



STD 



-q,0 



n (B4) 



with q being the number of the targeted quasi-particle 
subspace. Considering also the inverse product DT, the 
commutator [T, D] does not need to be calculated while 
evaluating dt,H^ m ' if 



m > max (o max ,T,D, O max>DT 



(B5) 



holds. 

As an example, we consider the second order calcula- 
tion (n — 2) given in Tabs. [T] and for the ground state 
energy (q = 0). Calculating diTT 2 ^ , the commutator of 
the monomials 



T 

D : 



i f J t t 



(B6a) 
(B6b) 



with O min (T) — O min (D) = 1 occurs. The numbers of 
local creation and annihilation operators are given by 



c T = 2 
c D = 2 



aT = 
a D = 2. 



(B7a) 
(B7b) 



In the normal-ordering of TD, no local operator can can- 
cel {std — 0) implying ctd = 4 and aTD — 2. For the 
product DT, sdt — 2 pairs of local operators may cancel 
implying ctd > 2 and otd > 0. Using Eq. (B4), we find 



maxJD = n - [2] - [1] = -1 
max,DT = U - \l\ - [0] = +1. 



(B8a) 
(B8b) 



Since the commutator [T, D] yields monomials with a 
maximum order of at most 1 in the calculation of dgH( 2 \ 
it can not yield relevant contributions. Hence it does not 
need to be evaluated at all. 



But in a calculation of order n > 2 or aiming at a 
higher quasi-particle subspace q > 0, Eq. (B4) yields 



higher upper bounds for the maximum order and thus 
the commutator must be evaluated explicitly. Note that 
this basic a-priori rule is only sensitive to changes of 
the quasi-particle numbers. It can not anticipate that 
the commutator in this example actually vanishes due to 
other properties of the hardcore algebra of the triplons. 



Appendix C: Extended A-Priori Simplification Rule 

The real-space structure of the commutator arguments 
T and D allows us to extend the above a-priori rule in 
analogy to the extended a-posteriori rule in App. [A] Let 
Ct and Cd be the clusters of the creation operators in T 
and in D, respectively. Analogously, At, Ad are the clus- 
ters of their respective annihilation operators. Normal- 
ordering the product TD can cancel local operators only 
on the intersection 

Std — At D Cd- (CI) 

Due to the locality of the triplon algebra, the commutator 
vanishes if none of the clusters overlap 



S 



TD 



'AS 



DT 



(C2) 



Thus the normal-ordered product TD definitely has local 
creation operators on the union cluster 

Ctd 2 C T U (C D \ S TD ) (C3a) 
and local annihilation operators on the union cluster 



A TD ^A D U (A t \ Std) 



(C3b) 



There may be additional creation or annihilation opera- 
tors, but no general statements can be made on their ex- 



istence. In this sense, the right hand sides of Eqs. ( C3a 



and ( C3b ) are minimum clusters for the normal-ordered 
product TD. They can be used in Eq. (A3) to obtain 



an upper bound for the maximum order O ma x,T_D and 
the corresponding reasoning is used to obtain O maXj DT- 
This makes it possible to avoid the computation of the 
commutator [T, D] . 

Moreover, one can use the intersections Std and Sdt 
to exploit the hardcore property of the triplons: The 
normal-ordered product TD will vanish if Ct and Cd \ 
Std are not disjoint or likewise if Ad and (At \ Std) 
are not disjoint because the creation or annihilation of 
two triplons is attempted on the same site. 

Although it is less strict, the basic a-priori has the ad- 
vantage to be much more lightweight in comparison to the 
extended a-priori rule because it requires mere counting 
of operators. Furthermore, it can be used very efficiently 
in the context of translation symmetry Because it does 
not rely on the real-space structure of a term, it can be 
applied to all terms in the translation group in contrast 
to the extended rule. Therefore, for best performance it 
turns out to be most efficient to combine both rules in 
practice. 
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Appendix D: Series expansion of ground state energy 
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Table IV. Coefficients of the perturbative evaluation for the ground state energy per rung for various y = 
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Appendix E: Series expansion of hopping terms 
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-0.99486640 


-0.27885442 


-0.35491362 


-0.33963371 


-0.11778655 


14 


-1.98783536 






1.36473269 




15 


-0.70791841 






2.41431850 






U 


h 


4 


-0.03906250 


-0.03116138 


-0.02792639 






5 


-0.04687500 


-0.03317747 


-0.03121078 


0.02734375 


0.01880795 


6 


0.03564453 


0.02548558 


0.01949955 


0.03906250 


0.02443938 


7 


0.13452148 


0.08115556 


0.07210783 


-0.03021240 


-0.01723325 


8 


0.08452606 


0.04273061 


0.04437608 


-0.13285828 


-0.06882735 


9 


-0.16891479 


-0.08803878 


-0.07128815 


-0.09355831 


-0.04418694 


10 


-0.37874413 


-0.17065003 


-0.15393594 


0.18216133 


0.07772144 


11 


-0.12521664 


-0.04650305 


-0.05271784 


0.44106736 


0.17211308 


12 


0.63520241 


0.24143238 


0.20728113 


0.16840086 


0.06173547 


13 


1.09758909 


0.37362081 


0.34061615 


-0.74875764 


-0.23909940 


14 


0.10661141 






-1.34708600 




15 


-2.19284584 






-0.12406771 






ts 


t 7 


6 


-0.02050781 


-0.01398812 


-0.01184453 






7 


-0.03417969 


-0.02091311 


-0.01819031 


0.01611328 


0.00921603 


8 


0.02478790 


0.01521149 


0.01095480 


0.03076172 


0.01599404 


9 


0.13085937 


0.06817047 


0.05621627 


-0.02009487 


-0.00957393 


10 


0.10456268 


0.04648122 


0.04328735 


-0.12947965 


-0.05581397 


11 


-0.19036049 


-0.08638434 


-0.06396735 


-0.11682585 


-0.04574993 


12 


-0.50204569 


-0.19572116 


-0.16328447 


0.19740278 


0.07034140 


13 


-0.20806881 


-0.06412749 


-0.07130687 


0.57210467 


0.18566781 


14 


0.90379496 






0.26874640 




15 


1.68652741 






-1.06182222 






ts 


tg 


8 


-0.01309204 


-0.00755457 


-0.00618286 






9 


-0.02819824 


-0.01468727 


-0.01219007 


0.01091003 


0.00519218 


10 


0.01595318 


0.00838267 


0.00574076 


0.02618408 


0.01132456 


11 


0.12834901 


0.05682685 


0.04493297 


-0.01227187 


-0.00490073 


12 


0.13062643 


0.05022786 


0.04326984 


-0.12730413 


-0.04569697 


13 


-0.20082120 


-0.07779488 


-0.05499657 


-0.14531211 


-0.04729919 


14 


-0.64685442 






0.20129516 




15 


-0.34698491 






0.72601478 






tio 


t 


i 


10 


-0.00927353 


-0.00450284 


-0.00360404 






11 


-0.02454758 


-0.01079578 


-0.00869977 


0.00800896 


0.00317279 


12 


0.00896406 


0.00409773 


0.00260778 


0.02318382 


0.00834536 


13 


0.12625622 


0.04714879 


0.03633844 


-0.00596249 


-0.00201914 


14 


0.16070462 






-0.12514982 




15 


-0.19870083 






-0.17663111 






ii2 


tl3 


12 


-0.00700784 


-0.00285462 


-0.00224969 






13 


-0.02202463 


-0.00814201 


-0.00643333 


0.00619924 


0.00204486 


14 


0.00321460 






0.02102351 




15 


0.12395082 






-0.00067934 






ti 4 


tis 


14 


-0.00553504 










15 


-0.02014753 






0.00498153 





Table V. Coefficients of the perturbative evaluation for the dispersion for y — 1 and y — 1.2. 
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